NON-LINEAR   ANALYSIS    OF    THE    SMOOTH    PURSUIT    EYE    MOVEMENT    SYSTEM    OF    THE    MACAQUE, 


by 
Lei and    S.^Stone    | 

Submitted  as  partial  requirement  for 
the  degree  of  Master's  of  Science 

at  the 

University  of  California, 
Berkeley,  California. 


Acknowledgements . 


The  author  would  like  to  thank  professor  Shankar  Sastry  for  his  Insight 
and  patience  both  of  which  were  essential  to  the  completion  of  this 
project.  He  is  also  indebted  to  professor  Steve  Llsberger,  Lauren  Westbrook, 
and  Ed  Morris  whose  experimental  results  were  graciously  made  available  to 
him  and  used  both  in  the  development  and  testing  of  the  model  presented 
here. 

• 
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ABSTRACT 


Using  the  time  domain  data  of  Llsberger  and  others,  a  nonlinear  model 
of  the  Smooth  Pursuit  system  of  the  macaque  was  constructed.  The  sinusoidal 
response  of  the  model  was  determined  using  a  modified  form  of  the  describing 
function  method.  The  results  were  compared  to  actual  frequency  domain  data 
in  one  monkey.  The  discrepancy  between  the  model's  response  and  the  actual 
data  is  Interpreted  as  being  the  result  of  a  prediction  process  which  is 
assisting  the  Smooth  Pursuit  during  predictable  tracking.   In  this  way,  a 
quantitative  estimate  of  the  role  of  the  predictor  is  made  at  least  for 
frequencies  between  .7  and  2  Hz  and  for  amplitudes  between  .3  and  2  deg. 


I.  Introduction  and  Background. 

In  the  macaque  and  all  foveate  animals,  the  -function  o-f  the 
oculomotor  system  is  to  put  the  image  o-f  an  object  o-f  interest  on 
the  -fovea,  the  region  o-f  the  retina  where  the  highest  resolution 
of  visual  processing  is  passible,  and  maintain  that  image  on  the 
•fovea  by  compensating  -for  any  movements  o-f  the  objects  or  the 
head  with  the  appropriate  eye  movements.  This  task  is  achieved 
by  several  subsystems  each  producing  a  speci-fic  category  o-f  eye 
movement:  the  saccadi_c_s^stem  which  produces  rapid  reor  i  entat  i  on 
o-f  gaze  is  used  to  scan  the  visual  environment  and  orrect  -for 
di -f -f  erences  between  the  object  (desired)  position  and  eye 

position,   the   vesti.bul^o-ocuLa.C reflex   produces    smooth    eye 

movements  equal  and  opposite  to  that  o-f  the  head  thereby 
stabilizing  images  during  head  turns,  and  the  smgpth_gursuit  (SP) 
system  produces  smooth  eye  movements  appropriate  -for  stabilizing 
the  image  o-f  a  moving  target  on  the  -fovea. 

To  characterize  what  aspects  of  the  visual  experience  are  used 
by  the  SP  system,  Rashbass  Cl]  per-formed  experiments  in  which 
position  and  velocity  error  were  in  con-flict  and  determined  that 
target  motion  not  target  position  is  the  relevant  input  to  the 
system.  It  was  concluded  that  the  SP  system  can  be  modelled  as 
simple  velocity  controller.  (Recent  evidence  has  shown  that 
Rashbass'  results  are  not  universally  applicable  and  that  in 
certain  situations  position  errors  can  also  produce  smooth  eye 
movements  [23  C3D  but  these  e-f-fects  are  small  and  not  relevant  in 
the  paradigms  used  below  so  will  be  ignored  here). 

On  purely  theoretical  and  logical  grounds  an  eye  velocity 
positive  -feedback  loop  was  postulated  C4D  to  account  -for  the 
accurate  steady-state  tracking  when  the  -feedback  velocity  error 
signal  is  small.  This  hypothesis  was  con-firmed  by  two 
experimental  results:  1)  in  both  monkey  and  man,  retinal  error 
velocity  is  linearly  correlated  with  eye  acceleration  during 
steady-state  sinusoidal  tracking  C5D,  and  2)  in  open  loop 
conditions,  constant  eye  velocity  is  maintained  in  the  absence  o-f 
visual  input  C63. 

The  SP  system  there-fore  can  be  considered  as  a  -feedback  system 
driven  by  image  motion  velocity  on  the  retina  and  whose  output  is 
eye  acceleration  used  to  modulate  ongoing  eye  velocity. 

Despite  the  1OO  msec  delay  in  the  system  response  C73,  this 
model  predicts  relatively  accurate  steady-state  tracking  of 
approximately  constant  velocity  targets:  the  SP  system  will 
accelerate  the  eye  up  to  the  appropriate  velocity  (the  saccadic 
system  will  then  correct  any  remaining  position  error)  and  hold 
it  there  by  modulating  eye  velocity  only  when  image  slippage 
occurs.  It  does  not  however  account  -for  the  accurate  tracking  o-f 
sinusoidal  motion  where  target  velocity  is  constantly  changing 


(predictably)  and  the  delay  should  but  does  not  cause  significant 
phase  lag  between  target  and  eye  velocity  C53.  A  predictor 
(target  trajectory  estimator)  is  there-fore  postulated  to  account 
•for  the  data.  Many  models  described  how  the  SP  system  could 
predict  C81  ,  but  as  yet  no  study  has  attempted  to  quantify  the 
actual  prediction  process  going  on  during  steady-state  tracking 
o-f  sinusoids. 

The  goal  o-f  the  present  study  is  to  use  in-formation  about  the 
open  loop  time  domain  response  of  the  SP  system  to  contruct  a 
model  that  will  be  used  to  estimate  the  theoretical  response  to 
sinusoidal  inputs  without  the  benefit  of  prediction.  The  model's 
response  can  then  be  compared  with  the  actual  open  loop 
sinusoidal  response  of  the  system.  In  this  way  a  quantitative 
measure  of  the  predictor's  contribution  to  the  frequency  response 
can  made  (at  least  within  the  frequency  and  amplitude  range 
examined).  It  is  hoped  that  the  use  of  modelling  to  determine 
the  overall  response  of  the  predictor  will,  in  combination  with 
future  experimental  results,  help  elucidate  how  the  SP  system 
overcomes  its  delay  and  dynamics  and  is  able  to  track  predictable 
targets  accurately. 


I  I.  The  Model . 


The  construction  of  the  model  was  based  on  the  results  o-f 
Lisberger  and  Westbrook  on  the  initiation  o-f  Smooth  Pursuit  eye 
movements  to  steps  in  target  velocity  C63.  They  found  that  the  SP 
system  responds  to  steps  in  target  velocity  by  producing  an  eye 
acceleration  with  a  latency  of  about  1OO  msec.  whether  the  eye 
was  initially  at  rest  < i ni ti ati on) or  v*s  tracking  a  constant 
velocity  target  (maintenance).  They  concluded  1)  that  the  1OO 
msec.  delay  in  the  tracking  response  represents  a  true  input 
time  delay  rather  than  a  "start-up"  delay  and  2)  that  the  initial 
1<~>0  msec  after  the  onset  of  eye  movement  can  be  considered  open 
loop  because  visual  feedback  is  also  delayed  by  1OO  msec.  <The 
later  hypothesis  was  tested  using  a  true  open  loop  paradigm  and 
found  to  be  correct.)  This  principle,  illustrated  in  the 
following  figure  (fig.  1),  was  used  to  study  the  open  loop 
performance  of  the  SP  system  as  a  function  of  several  parameters 
of  image  motion  (position,  velocity,  and  acceleration). 
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Fig.l  P>.  Schematic  illustrating  the  open  loop  interval  ?  B.  True 
open  loop  response  to  a  velocity  step  (stimulus  sho<nn>  with  the 
closed  loop  response  superimposed. 

They  first  studied  the  effect  of  initial  target  position  on  the 
response  to  velocity  steps  and  found  that  the  initial  10O  msec, 
could  be  subdivided  into  an  "early"  (~O-3O  msec.,  largely 
independent  of  initial  position)  and  a  "late"  (~3O-1OO  msec., 
strongly  spatially  dependent)  response.  They  therefore  divided 
the  1OO  msec.  open  loop  interval  into  2O  msec.  bins  and  the 
0-2O  msec.  bin  was  taken  as  a  measure  of  the  early  response 
while  the  6O-8O  msec.  bin  was  assumed  to  be  characteristic  of 


the  late  response. 

The  late  response  is  strongest  -for  targets  moving  near  and 
towards  the  fovea  (fig.  2A)  and  proportional  to  the  velocity  step 
size  (-fig.  2B)  up  to  around  45  deg/sec. 
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Fig. 2  f) .  Spatial  weighting  of  the 
B.  Magnitude  of  the  late  response 
size. 


late  response  (60- SO  msec  bin); 
as  a  function  of  velocity  step 


The  response  peaks  for  a  position  ®  between  0-3  deg.  off  the 
fovea  and  is  about  zero  for  positions  >1O  deg  if  the  target  is 
moving  away  from  the  fovea  and  for  positions  >2O  deg  if 
approaching  the  fovea.  This  spatial  weighting  of  the  late 
response  (the  output  of  element  A)  is  represented  in  the  model  by 
a  non-linear  function  with  hysteresis  (fig.  3). 


Fig.  3  Non-linear  element  f)  illustrating  the  hysteretic  nature  of 
the  spatial  weighting  of  the  late  response.   B  is  the  peak  angle 
parameter  C*0-3  deg). 


The  model  velocity  gain  (element  B)   is   constructed   from   the 


information  in  -figure  2B  and  -from  the  results  o-f  Morris  and 
Lisberger  C3DC93  that  show  that  the  low  velocity  (less  than  ~O-5 
deg/sec)  slope  is  much  higher  than  the  asymptotic  slope  (-fig. 
4A) .  (The  low  velocity  slope  m  varied  between  about  26.5  sec-1 
during  initiation  with  dim  background  illumination  and  33.33 
sec-1  during  maintenence  in  the  dark).  Figure  4B  illustrates  the 
non-linear  idealized  gain  element  used  in  the  model  -for  the  late 
or  velocity  pathway. 
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Fig. 4  ft.  Lou  velocity  responses;  B.  Non-linear  element  B  nith  ION 
velocity  slope  parameter  m  C*26-33  deg/sec/sec>  and  velocity 
elboN  parameter  fr  C"mrS  deg). 

The  early  response  is  only  weakly  spatially  weighted  (not 
shown)  but  is  strongly  related  to  the  acceleration  o-f  the  target 
at  the  onset  o-f  the  velocity  "steps"  (actually  velocity 
traperoids).  It  is  approximately  proportional  to  the 
acceleration  used  to  achieve  final  velocity  up  to  around  125 
deg/sec2  (-fig.  6A)  and  is  weighted  by  the  -final  velocity  up  to  15 
deg/sec  (-fig.  5A)  .  The  model  idealizations  o-f  the  gain  and 
velocity  weighting  o-f  the  early  or  acceleration  pathway  are  shown 
in  -figures  6B  and  5B  respectively. 
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5  ft.  Velocity  Weighting  of  the  early  response  (0-20  msec 
bin);  B.  Non-linear  element  C. 

Fig.  6  ft.  Early  response  (JO  msec  after  onset)  as  a  function  of 
the  acceleration  used  to  achieve  5  deg/sec  velocity;  B.  Non-linear 
element  D. 

In  -figure  7,  the  O  deg/sec  curve  represents  the  pursuit 
response  to  a  target  which  is  already  moving  at  constant  velocity 
when  pursuit  is  initiated  (pure  velocity  step)  while  the  in-f 
deg/sec  curve  represents  the  response  to  a  target  which  is 
initially  stationary  and  "instantaneously"  accelerates  up  to 
speed  at  time  zero  (velocity  step  plus  acceleration  impulse  at 
time  zero).  Assuming  dynamic  linearity  and  independence  (see 
below),  the  di-f-ference  between  the  two  curves  can  there-fore  be 
considered  as  the  acceleration  impulse  response  (fig.  7B)  and  the 
O  deg/sec  curve  can  be  considered  as  the  velocity  step  response 
(fig.  7A) .  The  dynamics  of  the  late  or  velocity  pathway  was 
determined  to  have  the  transfer  function  l/(.O4s  +  1)  and  that  of 
the  early  or  acceleration  the  transfer  function  l/(.Ols  +  1)". 
(The  responses  of  the  idealized  model  transfer  functions  are 
graphed  along  with  the  data  in  fig.  7). 
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Fig.  7  ft.  Velocity  step  rsponse  with  <  inf  deg/sec/sec>  and 
Mithout  (O  deg/sec/sec)  acceleration  impulse?  B.  Acceleration 
impulse  response. 

The  complete  model  is  given  in  -figure  8  with  parameters  8,  m,2f 
that  vary  within  the  constraints  given  above.  It  is  constructed 
such  that  the  input,  target  motion  as  a  -function  o-f  time  (where 
the  position  and  -first  two  derivatives  o-f  motion  are  the  only 
relevant  aspects  o-f  the  stimuli),  produces  an  eye  acceleration 
through  two  parallel  pathways.  The  model  assumes  that  the  time 
delay  is  at  the  input  end,  that  the  static  non-linearities 
discussed  above  -follow  with  the  weighting  parameters  (value  O-l) 
being  multiplicative  at  the  gain  output,  and  that  the  dynamics  o-f 
both  pathways  are  linear  (i.e.  independent  o-f  amplitude), 
independent  (i.e.  that  the  dynamics  are  independent  o-f  the  value 
o-f  the  weighting  parameter),  and  that  they  act  as  low  pass 
•filters  at  the  output  end  such  that  higher  harmonics  are 
removed.  (The  -first  two  assumptions  are  critical  to  the 
interpretation  o-f  the  data  in  -figure  7). 
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III.  Analysis. 

Open  loop  steady-state  sinusoidal  data  [10]  show  that  despite 
the  SP  system's  non-linearities  that  the  sinusoidal  response  for 
amplitudes  less  than  2  deg  and  -frequencies  between  .7  and  2  Hz.  is 
quasi  linear  i.e.  that  the  system  phase  shifts  and  scales  the 
sinusoidal  input  but  does  not  produce  much  harmonic  distortion. 
In  the  case  o-f  a  linear  system  the  gain  and  phase  are  -functions 
o-f  -frequency  only  whereas  -for  a  quasi  linear  system  they  may  be 
•functions  o-f  amplitude  as  well. 

This  -fact  suggested  that  the  method  o-f  describing  -functions 
(DF's)  might  be  appropriate  -for  the  sinusoidal  analysis  o-f  the  SP 
model.  This  method  merely  assumes  that  the  steady-state  periodic 
response  o-f  any  non-linear  element  written  as  a  Fourier  series 
can  be  truncated  after  the  first  term  without  significant  loss  of 
signal.  The  element  therefore  can  be  characterized  by  the 
amplitude  and  phase  of  the  first  harmonic  of  the  output  as  a 
function  of  input  amplitude  and  frequency.  The  DF  itself  is  a 
complex  gain  calculated  as  the  ratio  of  the  first  complex 
coefficient  in  the  Fourier  expansion  of  the  response  to  Asinwt 
divided  by  A.  In  the  case  of  static  nonl  i  near  i  ti  es,  it  is 
independent  of  frequency  and  given  by  the  following  integral 
CUD: 

V' 
0          '  (1) 


The  sinusoidal  analysis  of  the  SP  model  must  deviate  from 
classical  DF  theory  because  the  weighting  factor  non-linearities 
are  even  functions  and  therefore  have  no  first  harmonic  component 
in  their  output  (see  Appendix  A).  In  addition  their  effect  is 
multiplicative  so  that  not  only  the  D.C.  term  but  also  higher 
harmonics  must  be  considered  because  the  product  will  contain 
first  harmonic  terms  produced  by  the  multiplication  of  higher 
order  harmonics.  Specifically,  the  weighting  factor  elements 
produce  a  D.C.  term  and  even  harmonics  (aO,  a2,  a4,  a6.  .  .  )  whereas 
the  gain  elements  (odd  functions)  produce  no  D.C.  term  and  odd 
harmonics  (bl  ,  b3,  b5.  .  .  >  so  the  product  will  contain  first 
harmonic  terms  from  the  factors  aO*bl,  a2*bl,  a2*b3,  a4*b3, 
a4*b5,  etc.  In  this  case,  only  the  first  two  terms  of  this 
expansion  are  kept  because  the  higher  order  terms  are 
considerably  smaller  and  because  in  this  way  the  gain  elements 
could  be  analyzed  in  the  usual  DF  manner.  For  the  weighting 
factor  non-linearities  however  both  the  D.C.  gain  NO  and  second 
harmonic  complex  gain  N2  were  calculated  using  the  following 
equations  (derived  by  dividing  the  appropriate  Fourier 
coefficient  by  the  input  amplitude): 


IV.  Computation. 

(Details  of  the  computation  of  the  model  output  can  be  -found  in 
Appendices  A,  B  and  C) . 

The  classical  "first  order"  DF's  Mere  calculated  -for  the  odd 
non-linear  elements  B  and  D  using  the  integral  -formula  given 
above  (EQ.  1)  yielding! 


<  II  <  tf   <4> 

<  I  1  <  45 


<  12  <  125  (5) 


125  <  12 

with  II  the  amplitude  o-f  the  target  velocity  sinusoid  and  12  that 
o-f  the  acceleration  (see  -fig.  8).   I-f  A  and  -f  are  the  amplitude  and 
•frequency  o-f  the  actual  input  target  position  sinusoid  then 
Il=2*PI*f*A  and  12  =  2  *  PI  *  -f  *  1  1  .    and  m  are  de-fined 
in  -fig.  4. 

The  -first  order  DF's  -for  the  even  nonl  i  near  i  ti  es  A  and  C  are 
zero.  The  zeroth  order  (D.C.  )  and  second  order  terms  were 
calculated  using  equations  (2)  and  (3)  respectively  (see  Appendix 
A)  . 


0  <  A 


<  A  <  10 


defined  in 


The  Bode  Plot  of  interest  is  that  of  eye  velocity  E  with 
respect  to  target  velocity  T  because  of  the  available  data  for 
verification  of  the  model  (see  next  section)  therefore  the  phase 
of  t  will  be  defined  as  zero  and  all  phases  will  be  calculated 
with  respect  to  this. 


The  n-th  order  (n  >  O)  output  of  the  non-linear  elements  were 
calculated  using  the  -following  formula  where  np  and  nq  are  the 
real  and  imaginary  parts  o-f  the  n-th  order  DF  respectively: 


Input  Amp.  *1r\6f  + A<»X  #sin(nwt  -  arctan  (np/nq)  •+•  (PI/2)  +  input  phase)  (B) 
which  when  the  DF  is  real  valued  simplifies  to: 


Input  Amp.  *  np  *  sin(nwt  +  input  phase) 

and  the  D.C.  term  (zeroth  order  output)  is  given  by: 

Input  Amp.  *  Zeroth  Order  DF 

The  above  formulae  were  used  to  calculate  the  outputs  of  the 
various  non-linear  elements:  YO,  Y2  «*>  ,  Xl(f),  WO,  W2  (/3  >  ,  X2(^>. 
(The  notation  >:  (y)  meaning  ampl  itude  (phase)  and  the  variables 
defined  in  fig. 8).  The  outputs  of  the  multiplier  stages  were 
calculated  using  standard  trigonometric  identities  (see  Appendix 
C).  The  D.C.  term  of  the  weighting  factor  multiplies  the  output 
of  the  gain  element  yielding  a  first  order  term  and  the  second 
order  term  of  the  weighting  factor  also  produces  a  first  order 
term.  The  sum  of  these  two  first  order  output  terms  is  the 
output  of  the  multiplier:  Viith 

'o^  ^x(^)  N*  OC  4  X 

vi  ttr«*V 

c.s  ».  *  . 


(9) 


(10) 


The  two  multiplier  output  terms,  Z1(.)  and  Z2(ty  are  then 
filtered  linearly  yielding  01  (J,)  and  O2  <£)  ,  the  outputs  of  the 
velocity  and  acceleration  pathways  respectively.  These  two  are 
then  added  yielding  the  output  eye  acceleration  which  is 
integrated  to  get  eye  velocity. 


(11) 


The   calculation   of  the  model  output  was  done  on 
using  the  program  Ismodel  written  in    programming 
copy  of  the  program  is  in  Appendix  B. 


a   PDP-11/23 
language  C.  A 


V.  Results  and  Discussion. 


The  results  of  the  simulation  are  shown  in  Table  1.  The 
parameters  8,  m,  and  f ,  were  optimired  by  eye  within  the 
constraints  set  -forth  in  section  II.  The  third  and  -fourth  columns 
are  the  output  of  the  model  in  -fig. 8,  the  -fifth  and  sixth  columns 
the  output  o-f  the  model  i -f  the  early  or  acceleration  pathway  is 
ignored,  and  the  seventh  and  eighth  columns  the  experimental 
results  o-f  Lisberger  (1O)  o-f  the  actual  open  loop  -frequency 
domain  response  o-f  one  monkey.  The  optimization  procedure  was 
not  simply  a  curve  -fitting  procedure  because  1)  the  allowed 
values  were  severely  constrained  -from  the  experimental  data,  2) 
the  number  o-f  parameters  (3)  did  not  offer  enough  degrees  of 
freedom  for  any  real  "manipulation"  of  the  output,  3)  the 
simulation  was  insensitive  to  any  change  in  the  value  of  0  within 
its  constrained  value  interval,  and  4)  the  optimal  values  for  m 
and  within  the  constrained  bounds  were  very  close  to  the  best  a 
Br.i9r.i  estimates.  It  was  performed  more  to  determine  the  effect 
of  the  parameter  modulation  on  the  performance  than  to  match  the 
data. 

As  mentioned  above,  the  model  performance  is  insensitive  to 
changes  in  0  so  it  is  difficult  to  choose  a  value  with  any 
conviction  but  the  data  of  Morris  (9)  suggests  that  it  may  be  as 
low  as  .5  deg  or  lower  so  the  value  of  .5  was  used  in  the 
simulation  in  Table  1.  The  low  velocity  slope  m  of  the  velocity 
gain  element  has  a  value  that  varies  between  26.5  /sec  for 
initiation  of  pursuit  with  dim  background  and  33.3  /sec  for 
maintenance  (modulation  of  an  ongoing  eye  velocity)  of  pursuit  in 
the  dark  C33C9D.  The  a  priori  best  estimate  is  therefore  33.3 
given  that  the  conditions  under  which  it  was  derived  most  closely 
match  those  under  which  the  frequency  domain  data  was  found.  The 
optimal  value  was  found  to  be  32  /sec.  Increases  in  m  had  the 
effect  of  almost  uniformly  raising  the  gains  for  all  amplitudes 
and  frequencies.  The  velocity  gain  curve  bends  starting  at 
around  1.5  deg/sec  and  achieves  a  new  constant  slope  at  around  5 
deg/sec.  Most  of  the  gain  attenuation  occurs  at  the  low  end  of 
the  interval  so  the  best  a  priori  guess  would  be  somewhere 
between  1.5  and  3  deg/sec.  The  optimal  value  of  "o  was  determined 
to  be  2.75  deg/sec.  Increasing  If  had  the  effect  of  decreasing 
the  high  frequency  and  amplitude  attenuation  in  gain.  One 
encouraging  result  of  this  is  that  the  optima  found  within  the 
constrained  intervals  seemed  also  to  be  the  global  optima  in  that 
values  outside  the  intervals  may  the  simulation  less  accurate. 

As  can  be  seen  from  the  results  in  Table  1  and  fig.  9,  the 
model  relatively  accurately  predicts  the  gain  data  exept  that  the 
high  frequency  attenuation  in  the  model  occurs  at  somewhat  too 
low  a  frequency  (at  <1  rather  than  >1  Hz)  and  is  somewhat  too 
broad  in  nature.  No  manipulation  of  the  parameters  could  improve 
the  overall  shape  of  the  gain  curve.  The  phase  data,  as 
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expected,  is  not  at  all  predicted  by  the  model.  One  very 
interesting  result  however  is  that  the  slope  of  the  phase  curve 
is  approximately  correct  especially  -for  the  higher  amplitude 
curves  and  that  there  is  a  constant  phase  discrepancy  o-f  about 
140  deg  -for  all  frequencies  and  amplitudes  examined. 

It  should  also  be  noted  that  the  acceleration  pathway 
contributes  very  little  to  the  response  as  can  be  seen  -from  a 
comparison  o-f  columns  -four  and  six,  and  -five  and  seven.  In  an 
e-f-fort  to  examine  the  potential  role  of  the  acceleration  pathway 
(assuming  the  gain  estimate  used  in  the  model  construction  was 
possibly  inaccurate  and  too  low),  the  acceleration  gain  and 
saturation  point  were  manipulated  arbitrarily  while  attempting  to 
improve  the  fitting  of  the  data.  <A  modified  form  of  the  program 
Ismodl  was  used  where  element  D  of  the  model  was  modified).  The 
results  are  shown  in  Table  2.  When  the  acceleration  gain  is 
increased  by  a  factor  of  ten  and  the  saturation  acceleration  by  a 
factor  of  two,  the  phase  lag  can  be  reduced  to  the  point  where 
the  1.66  deg  2  Hz  model  response  matches  the  data  almost  exactly 
but  the  shape  of  both  the  gain  and  phase  curves  is  seriously 
compromised.  No  systematic  optimization  was  performed  so  it.  is 
difficult  to  make  any  firm  conclusions  particularly  since  the 
parameter  manipulation  for  the  acceleration  pathway  was 
arbitrary.  It  can  however  be  concluded  that  the  pathway  can  be 
made  to  contribute  substantial  phase  lead  to  the  response  and 
therefore  may  potentially  be  involved  in  the  prediction  process. 

If  the  premise  that  the  model  accurately  represents  the  SP 
system  sans  prediction  were  valid,  the  finding  that  the  model 
lags  the  data  by  around  14O  deg.  at  all  frequencies  and 
amplitudes  is  an  interesting  one.  Any  future  model  of  the 
predictor  must  account  for  this  i.e.  the  predictor  adds  140  deg 
worth  of  phase  lead  to  the  system  with  little  effect  on  the  gain 
over  the  amplitude  and  frequency  range  examined. 

In  the  future,  other  model  validation  experiments  would  be 
beneficial:  1)  model  simulation  can  be  used  for  velocity  step 
inputs  at  different  positions  in  an  effort  to  determine  if  the 
model  can  reproduce  the  data  that  was  used  to  derive  it  (i.e.  is 
the  model  faithful?)  and  2)  as  has  been  pointed  out  in  C123,  the 
ordering  of  the  non-linear  elements  and  linear  systems  in  each 
pathway  of  the  model  is  important  (i.e.  non-linear  elements  and 
linear  systems  do  not  commute)  and  multitone  experiments  could  be 
used  to  determine  if  the  proposed  ordering  in  the  model  is  a  good 
one.  These  experiments  would  provide  an  independent  method  of 
verifying  the  modelling  procedure.  In  addition  ,  if  the  model 
could  successfully  simulate  the  response  to  a  variety  of 
nonpredi cti ve  target  trajectories  then  the  premise  used  above  to 
determine  the  role  of  the  predictor  in  sinusoidal  tracking  would 
be  proven. 

Once  one  is  confident  of  the  model's  validity,  one  could  use  it 


in  combination  with  hypothesized  prediction  mechanisms  (using 
either  visual  and/or  occulomotor  information  to  construct  an 
estimate  o-f  target  trajectory  T  that  is  used  as  the  input  to  the 
model)  in  an  attempt  to  match  the  sinusoidal  data.  Experimental 
studies  o-f  the  open  loop  responses  to  predictable  targets  would 
be  used  in  an  e-f-fort  to  uncover  the  structure  o-f  the  prediction 
process.  These  results  would  then  be  incorporated  into  the 
predictor  model.  It  is  hoped  that  in  this  way,  a  predictor  model 
could  be  -found  that  accurately  mimicked  the  data  while  using  a 
physiologically  plausible  mechanism. 
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ILJU  5  B-l 

iclude  <std.h> 
iclude  <rtll.h> 

static  double  aC63  =  { . 32, . 42, . 62, . 83, 1 . 25, 1 . 66>  ; 

static  double  fC33  ={.7,1.,2.>; 

static  double  gdataC63C33  =  C6.  13,  5.  32,  .  96,  0.  ,  4.  9,  .  59,  0.  ,  6.  21 ,  1 .  29, 

5.34,5.35, 1.4,O. ,4.61, 1.36,3.87,3.71, 1. 17> ; 
static  double  pdata[63C33  «O7.  1 , -47.  7, -38.  5,  0.  , -36.  7, -45.  6,  0.  , -19.  3, 

-42.2,22.4,6. ,-55.2,O. , -16. 4, -4O. 5, 35. 1, 14.5, -49.  33- ; 
double  il  0,  i2  0,  wO  0,  w2  O,  xl  0; 
double  >:2  0,  yO  O,  y2  0,  zl  O,  z2  0,  ol  0,  o2  0; 
double  alpha  0,  beta  0,  theta  0,  eta  0  ,  slope  0; 
double  phil  O,  phi2  O,  psi 1  O,  psi2  0,  ki 1  O,  ki2  Oj 
double  nap  0,  naq  0; 
int  i  O,  j  0; 

double  gainC63C33  C0> ,  phaseC63C33  CO,  gC63C33  <03  ,  pC63C33  <0>; 
double  sqrtO,  cos  ( )  ,  sin(),  asinO,  atanO,  atan2(),  power  ()  ; 

Lain  ( ) 

put-fmt  ("Peak    Angle    Parameter     (theta)     ?    "); 

putch(-O12) ; 

get-fmt  ("V.d",     ?<theta>  ; 

put-fmt  ("Velocity    Elbow    Parameter     (eta)     7    "); 

putch(-012) ; 

getfmt  ("*/.d",    ?<eta)  ; 

put-fmt  ("Low    Velocity    Slope    7    "); 

putch(-012) ; 

get-fmt  (  "7.d"  ,  Aslope)  ; 

•for  (  i  =  O  ;  i  <  6  ;  i++> 

•for  (  j  =  0  ;  j  <  3  ;  j++) 
< 

/*  Velocity  Pathway  */ 

il  =  2  *  3.1415926  *  fCj]  *  aCiD; 
phil  =  -  .2  *  3.1415926  *  fCjD; 

/*  Velocity  Gain  Element  */ 
i-f  (i  1  >  eta) 

xl  =  (.6366  *  (slope  -  6.2222) 

*  (asin(eta/il)  +  (eta/il)  * 

sqrt  ( 1  -power (eta/i 1,2.))) +6. 2222) *i 1 1 
el  se 

x 1  =  si  ope  *  i 1 ; 

/*  Spatial  Weighting  */ 

i-f  (  aCi  3  >  theta) 

{        yO  =  ( (2OO  -  25  *  theta)  *  3.1415926  + 

3O  *  theta  *  asin(theta  /  aCi:>  + 
30*aCi  3*sqrt (1  -power (theta/aCi  3,2. ) ) ) 
/  (3.1415926  *  (10  +  theta)  * 
(20  -  theta) ) : 
nap  =  (2  *  (-30  *  theta  *  (1  - 

power (theta/aCi 3,2. ))  - 

.66667  *  aCi3  *  (10  -  2  *  theta  +  3O  * 
power (theta/aCi3,3. ) ) ) )/(3. 1415926*aCi3 

*  (10  +  theta)  *  (2O  •  theta)); 

naq  =  (2O  *  sqrtd  -  power  ( theta/aC  i  3  ,  2.  )  )  * 
(1  +  5  *  power (theta/aCi 3,2. )))  / 
if  1  /n  =:O->A*  r  1  o  -  thDta)  t  (?ft  -  theta)): 


y^  =  sqrt  (power  (nap,  ^'.  )-«-power  (naq,^.  )>  *  aiij; 

alpha  =  phil  -  atan2 (nap, naq) ; 
3- 

else 
<        yO  =  10  /  (10  +  theta) ; 

y2  =  (8  t  aCiD)  /  (3  t  (10  +  theta) 
t  3. 1415926) ; 

alpha  =  phil  +  1.5708; 


/*  Crass-Product  Calculation  */ 

zl  =  xl  *  sqrt (power (yO, 2. )  *  . 25  *  power (y2, 2. ) +  yO  * 
y2  *  cos(2tphil  -  alpha    1.5708)); 

psil  =  1.57O8  -  atan2(2  *  yO  *  cos  (phil)  •»•  y2  * 

cos(alpha  -  phil  +  1.5708), 2  *  yO  *  sin(phil) 
+  y2  *  sin(alpha  -  phil  +  1.5708)); 

/t  Linear  Filtering  #/ 

ol  =  zl  /  sqrt(power(2  *  3. 1415926** C j 3 , 2. ) *. 0016  +1 ) ; 

kil  =  psil  -  atan(.08  *  3.1415926  *  -f  C  j  3 )  ; 


/t  Acceleration  Pathway  t/ 

i2  =  2  *  3.1415926  *  i C j 3  *  il; 
phi  2  =  phil  +  1.57O8; 


/*  Acceleration  Gain  Element  #/ 
i-f  (12  >  125) 

x2  =  1.28  *  .6366  * (asi n ( 125/i 2)  +  (125/i2)  * 

sqrtd.  -  power  (125./i2,2.  )  )  )  *  i2; 
else     x2  =  1.28  *  i2; 


/*  Velocity  Weighting  */ 

if  (il  >  15) 

{        wO  =  (3.1415926  -  2*asin ( 15/i 1 )  +  .1333  * 

il  *  (1  -  sqrtd  -  power  (15/i  1,2.  ))))/ 
(3. 1415926  til); 
w2  =  1.2732  *  ((.O222  -  5  /  power (i 1 , 2. ) )  * 

sqrtd  -  power  (15/i  1,2.  ))-  .O222)  til; 
i-f  (w2  >  0) 

beta  =  phi  1  +  1.5708; 
else     beta  =   phil  -  1.5708; 
> 

else 

{  wO  =  . O4244  til; 
w2  =  .02829  til; 
beta  =  phil  -  1.5708; 


/t  Cross-Product  Calculation  t/ 

z2  =  x2  t  sqrt (power (wO, 2. )  +  . 25  t  power (w2, 2.)  + 
wO  t  w2  t  cos(2  t  phi2  ••  beta  -  1.5708)); 

psi2  =  1.5708  -  atan2(2  t  wO  t  cos  (phi  2)  -»-  w2  t 
cos (beta  •  phi  2  *  1.5708) ,2  t  wO  t 
sin(phi2)  +  w2  t  sin(beta  -  phi2  *  1.57OB)); 


/t  Linear  Filtering  t/ 

o2  =  z2  /  power  ( .  OO3948  t  -fCj]  t  -FCj3  *  1,2.); 


/*  uutput  calculation  */  j« 

gCi  3C  j3  =  ol/i2;  /  "*  * 

pCiDCjD  =  57.29578  *  (kil    1.5708)} 

gainCi3Cj3  =  sqrt (power (ol , 2. ) +power (o2, 2. ) +  2*ol*o2  * 

cos(ki 1  -  ki2) )  /  12; 
phaseCiDCjD  =  57.29578  *  (    atan2(ol  *  cos(kil)  + 

o2  *  cos(ki2),ol  *  sin(kil)  +  o2  *sin(ki2)))| 


/*  Printout  */ 

eptext ("Model  Output  -for  Peak  Angle  Parameter  Value  Theta  =  "); 

epdoub (t beta, 3, 2) ; 

epchar  ( ' \n' ) ; 

eptext ("Velocity  Elbow  Parameter  Eta  =  "); 

epdoub (eta, 3, 2) ; 

epchar  C\n'  )  ; 

eptext ("Low  Velocity  Slope  =  "); 

epdoub (slope, 5,3); 

epchar  ( ' \n' ) ; 

epchar  ( ' \n' ) ; 

eptext ("Amp. ") ; 

epchar  C\t'  )  ; 

eptext ("Freq. ") ; 

epchar  C\t'  )  ; 

eptext ("Gain") ; 

epchar (' \t' ) ; 

eptext ("Phase") ; 

epchar  C\t'  )  ; 

eptext ("Gain") ; 

epchar <'\t' ) ; 

eptext ("Phase") ; 

epchar  C\t'  )  ; 

eptext ("Gain") ; 

epchar  C\t'  )  ; 

eptext ("Phase") ; 

epchar ( ' \n' ) ; 

epchar  (' \t' ) ; 

epchar  C\t'  )  ; 

eptext ("Model ") ; 

epchar  C\t'  )  ; 

eptext ("Model ") ; 

epchar  C\t'  )  ; 

eptext ("No  Ace. ") ; 

epchar  C\t'  )  ; 

eptext ("No  Ace. ") ; 

epchar  C\t'  )  ; 

eptext ("Data") ; 

epchar  C\t'  )  ; 

eptext ("Data") ; 

epchar  ( * Xn' ) ; 

epchar (' \n' ) ; 


•for  (i  =  0  ;  i  <  6  ;  i++) 

•for  (j  =  0  ;  j  <  3  ;  j++) 
<        epdoub (aCi 3,5,2) ; 

epchar  C\t'  )  ; 

epdoub (f [ j  3 , 3. 1 ) ; 

epchar (' \t" ) ; 

epdoub (gainCi  3  L  j  3 , 5, 2) ; 

epchar (' \t7 ) ; 

pnrlnuh  (ohaspC  i3Ci3.4.1) 


epdouib  (gLiJLjJib,^); 

epchar  <'\t'  )| 

epdoub  (pC  i  ]  L  j  3  ,  4  ,  1)  ; 

epchar  C\t'  )  ; 

i-f  (gdatati  DC  jD  '=  0.  ) 

{        epdoub  (gdataCi  DC  jD,  5,2)  ; 

epchar  C\t'  )  5 

epdoub  (pdataCi  DC  j  3,4,  1)  } 

epchar  C\n'  )  ; 


else 


epchar  C\n'  )  ; 


epchar  ( '  \-f '  )  ; 
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